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1.  Introduction 

The  specification  of  prognostic  variable  values  on  the 
lateral  boundaries  of  limited-area  numerical  models 
continues  to  be  a  serious  problem  in  atmospheric  modeling. 
The  computational  model  requires  boundary  values  where 
natural  atmospheric  boundaries  do  not  exist.  Theory  has 
provided  little  guidance  in  the  method  of  choosing  these 
boundary  values.  For  this  reason,  lateral  boundaries  have 
become  an  area  in  which  pragmatic  experimentation  has  been 
the  primary  source  of  boundary  condition  formulation.  This 
is  especially  true  in  light  of  theoretical  results  which 
show  that,  in  the  general  case,  it  is  impossible  to 
formulate  truly  well-posed  boundary  conditions  (Klemp  and 
Lilly,  1978).  A  large  number  of  lateral  boundary  conditions 
have  been  formulated  and  used  successfully  in  various  types 
of  atmospheric  models.  No  particular  formulation  has  been 
shown  to  be  clearly  superior,  and  indeed,  the  "best" 
boundary  condition  is  usually  determined  by  the  nature  of 
the  problem  under  consideration. 

The  problem  under  consideration  in  this  study  is  the 
successful  modeling  of  the  meso-beta  scale  with  a 
three-dimensional  hydrostatic  model.  The  lateral  boundary 
conditions  are  particularly  crucial  on  this  scale  because 
the  model's  grid  spacing  is  small  enough  to  resolve 


high  frequency  atmospheric  waves  (such  as  internal  gravity 
waves  amd  lee  waves)  while  the  time  and  space  scales  are 
large  enough  to  require  inclusion  of  the  Coriolis  effect. 

A  successful  boundary  condition  must  minimize  the  reflection 
of  outgoing  waves  by  the  boundary  since  the  reflected  wave 
energy  might  destroy  the  interior  solution.  The  boundary 
condition  must  also  accurately  prescribe  the  boundary 
values  or  else  the  errors  may  be  geostrophically 
tramsmitted  to  the  interior. 

A  complete  discussion  of  geostrophically  adjusted 
boundary  errors  was  presented  by  Anthes  and  Warner  (1978). 

We  can  think  of  the  mean  values  of  velocity  ajid  temperature 
over  the  domain  in  terms  of  wavenumber  zero  flow.  The 
mean  velocities  in  the  domain  are  related  geostrophically 
only  to  gradients  of  the  mass  variables  across  the  domain. 

An  error  in  one  of  the  mass  variables  on  the  boundary  will 
imply  a  change  in  the  large  scale  gradient  and  lead  to 
a  change  in  the  mean  geostrophic  velocity  in  the  domain. 

The  mean  geostrophic  velocity  error  is  proportional  to  the 
boundary  value  error  divided  by  the  domain  width.  Thus, 
the  meso-beta  scale  should  be  much  more  suseptable  to  this 
type  of  error. 

This  report  will  present  several  lateral  boundary 
condition  formulations  in  the  context  of  the  meso-beta 
scale.  The  emphasis  will  be  on  radiation-type  conditions 
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because  it  will  be  shown  that  this  type  holds  the  most 
promise  of  acceptcUale  results.  This  study  is  not  meant  to 
represent  an  exhaustive  comparison  of  boundary  condition 
formulations,  but  it  does  provide  a  comparison  of  a  large 
number  of  conditions  and  especially  considers  several 
variations  on  the  radiation  condition. 

For  computational  economy,  ease  of  interpretation,  and 
availcQjility  of  theoretical  results,  the  study  reported  here 
concentrates  on  the  use  of  a  one-dimensional,  shallow  water 
model  to  test  the  various  boundary  condition  formulations. 
This  model  filters  out  many  forms  of  waves  important  in  the 
atmosphere  and  cautions  about  the  generality  of  some  of  the 
results  obtained  will  be  made  at  the  end  of  this  report. 

The  model  provides  a  good  test-bed  for  the  boundary 
conditions,  however,  and  yields  some  very  important  results 
concerning  modeling  on  the  meso-beta  scale. 

2.  The  One-Dimensional  Model 

The  one-dimensional  equations  governing  the  velocities 
u  and  V  and  the  depth  (see  fig.  1)  using  the  shallow 
water  approximation  on  an  f -plane  can  be  written 

af(i7v)  +  |7(ui7v)  +  ,  ^2) 

d?  *  0  .  <  3 ) 
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where  g  is  the  acceleration  of  gravity,  f  is  the  Coriolis 


parameter,  and  v  is  an  eddy  viscosity.  The  height- 


weighted  velocities,  tju  and  17V,  are  analogous  to  the 


pressure-weighted  velocities  used  as  prognostic  variables  in 


many  meso-scale  hydrostatic  models  (Anthes  and  Warner,  1978; 


Nickerson,  1979)  and  the  continuity  equation,  (3),  is 


analogous  to  the  pressure  tendency  equation  in  a  hydrostatic 


model.  In  fact,  the  set  (1)  through  (3)  is  identical  to  the 


equation  set  which  is  obtained  when  a  hydrostatic  mesoscale 


model  is  collapsed  to  one-dimension  along  the  x-axis. 


We  nondimensionalize  the  equations  using  the 


undisturbed  height,  H,  as  a  scale.  Thus,  letting  primes 


denote  the  nondimens ional  quantities,  we  have 


77  =  H77'  , 
u  =  \/gH  u'  , 

V  =  V'  , 

t  =  H/>/gH  t'  t'  , 


X  =  Hx' , 


f  =  /gFT  /Hf  '  =  x/gTH"  f '  , 
V  =  HVgH  K ' . 


Then  the  nondimensional  equations  are  (dropping  primes) 


f7(77u)  +  |^(u77u)-  f(77v)  = --^1^  +  , 

d?  +  *  0 


Note  that  K  may  now  be  regarded  ae  an  inverse  Reynolds 
number . 

The  choice  of  a  value  for  the  nondimensional  Coriolis 
parameter  determines  a  relationship  between  the  height  scale 
and  the  implied  latitude.  Thus,  if  a  latitude  of  interest 
is  prescribed,  a  dimensional  length  scale  must  be  implicitly 
chosen  even  in  the  nondimensional  analysis.  We  chose  to  let 
f  =  10”^  3“'  (  <p  =  43®Lat),  and  choose  a  height  scale  H  =  8 
lun  (approximately  the  density  scale  height).  This  yeilds  a 
gravity  wave  speed  of  \/gH  =  280  m/s. 

The  equations  are  finite  differenced  on  a  staggered 
grid  as  shown  in  fig.  1.  Time  integration  is  accomplished 
with  the  leapfrog  scheme  with  and  Asselin  filter  (Asselin, 
1972)  to  control  timesplitting.  The  Asselin  filter  factor 
is  set  at  0.5  as  suggested  by  Schlesinger  et  al.  (1983)  for 
the  "lagged”,  weak  diffusion  used  in  the  model.  The  normal 
model  domain  has  26  grid  points  for  17.  A  larger  domain 
with  50  grid  points  for  17  can  also  be  run.  The 
nondimensional  grid  spacing  is  chosen  to  be  Ax  =  2.5  which 
corresponds  to  20  km.  This  value,  and  the  resulting  domain 
size  of  500  km  is  typical  of  meso-beta  models.  It  also 
assures  that  all  waves  resolved  in  the  model  satisfy  the 
shallow  water  assumption  on  which  the  equations  are  based. 
With  the  parameter  values  chosen,  a  nondimensional  time  of 
1.0  corresponds  to  28.6  s.  The  timestep  is  chosen  to  be 


At  =  0.8.  This  is  aihout  half  the  size  required  for  linear 
stability  (including  the  reduction  required  by  the  use  of 
the  Asselin  filter) ,  but  is  typical  of  the  value  required  in 
a  2-D  or  3-D  model  with  the  same  Ax  value.  The  inverse 
Reynolds  number  is  set  at  K  =  7.813  x  10“^  which 
corresponds  to  an  eddy  viscosity  of  1.75  x  10^  m^/s.  The 
model  has  been  coded  with  a  variety  of  lateral  boundary 
conditions.  These  will  be  described  in  the  following 
sections  along  with  the  results  of  each  type  of  condition. 

3.  Description  of  the  Adjustment  EIxperiments 

The  major  concern  of  most  lateral  boundary  condition 
formulations  is  their  ability  to  allow  waves  generated  in 
the  domain  to  exit  (or  be  ebsorbed)  without  reflection.  If 
this  is  accomplished,  the  lateral  boundary  will  appear 
transparent  and  allow  the  limited  area  model  to  better 
represent  the  unbounded  region  it  is  simulating.  Another 
concern  is  the  possibility  of  errors  in  the  large  scale 
gradients  implied  by  the  boundary  values  being 
geostrophically  adjusted  throughout  the  domain.  On  the 
meso-beta  scale,  temperatures  and  windspeeds  must  be 
specified  with  errors  of  less  than  about  1°C  and  2  m/s, 
respectively,  in  order  to  give  acceptable  geostrophically 
adjusted  interior  errors  (Anthes  and  Warner,  1978). 

A  major  source  of  waves  in  the  domain  during  a 


simulation  is  the  geostrophic  adjustment  of  imbalances 
between  the  dynamic  fields.  In  general,  at  the  start  of  the 
integration,  the  velocity  and  mass  fields  will  not  be  in 
perfect  gradient  balance  and  gravity  and  Lamb  waves  will  be 
generated  during  the  adjustment  process.  These  waves  radiate 
away  from  the  imbalance  to  leave  behind  a  balanced  state. 

If  the  lateral  boundary  conditions  of  a  numerical  model 
reflect  some  or  all  of  the  wave  energy  back  into  the 
interior,  the  reflected  waves  will  unrealistically  influence 
the  simulated  solution,  interfering  with  the  balance  and 
prohibiting  the  adjustment  process. 

In  this  study,  experiments  were  conducted  with  two 
different  types  of  initial  imbalances  as  initial  conditions. 
Initial  condition  1  is  a  perturbation  in  the  height  field  and 
initial  condition  2  is  an  imbalanced  local  wind.  The  scale 
of  these  imbalanced  regions  is  smaller  than  the  domain  size 
of  the  model  (500  icm)  and  therefore  much  smaller  than  the 
Rossby  radius  of  deformation  (X  =  C/f  =  2.8  x  10®  km). 

Theory  predicts  that  the  pressure  should  adjust  to  the  winds 
with  the  wind  field  remaining  essentially  unchanged. 

Each  of  these  two  initial  conditions  were  ru.i  in  the 
large  domain  model  in  order  to  determine  their  natural 
evolution  during  adjustment.  For  these  simulations,  the 
model  used  a  radiation  lateral  boundary  condition  with  the 
phasespeed  set  to  the  analytical  value  of  C  =  1.0.  This 


will  be  shown  to  be  the  least  reflective  formulation. 

The  first  imbalance  (initial  condition  1)  is  a 
perturbation  in  height  only  with  the  velocities  set  to  zero. 
The  initial  height  field  and  its  evolution  are  shown  in  fig. 
2.  The  perturbation  is  a  height  deficit  of  0.3  % 

(comparable  to  a  pressure  deficit  of  3  mb)  in  the  center  of 
the  domain.  Height  excesses  are  located  on  either  side  of 
the  deficit  so  that  the  equalibrium  height  is  17=  1.0.  As 
linear  shallow  water  theory  would  predict,  the  evolution  of 
this  disturbance  is  for  two  deficit  waves  to  propagate  in 
opposite  directions  at  the  shallow  water  wavespeed,  C  =  ±1.0 
leaving  a  (trivially)  balanced  state  behind.  The  Coriolis 
effect  plays  little  role  in  this  adjustment  as  indicated  by 
a  simulation  with  f  =  0  (not  shown)  which  yeilded 
essentially  identical  results.  The  primary  velocities  are 
in  the  u  component  due  to  the  height  gradient.  Small  values 
of  V  are  accelerated  geostrophically  as  the  first  half  of 
the  wave  passes  a  point.  However,  these  are  quickly 
diminished  as  the  sign  of  the  acceleration  changes  during 
the  passage  of  the  rest  of  the  wave.  Note  that  by  t  =  72, 
the  wave  is  no  longer  present  in  the  area  covered  by  the 
small  domain  and  the  height  field  is  level. 

In  the  second  adjustment  experiment,  the  height  field 
is  initially  specified  as  rj  -  1.0,  but  a  limited  jet  is 
imposed  in  v  with  u  initially  zero.  This  is  essentially  the 


geo3trophlc  adjustment  problem  studied  by  Cahn  (1945)  using 
a  linear  shallow  water  model.  The  t  =  0  frame  of  fig.  3 
shows  the  extent  of  the  jet,  which  was  set  at  v  =  10  m/s 
(dimensionally).  The  remainder  of  the  figure  shows  the 
evolution  of  the  height  field.  As  predicted  by  linear 
theory,  the  height  field  adjusts  to  the  imbalance  by 
producing  a  wave  of  excess  which  moves  to  the  right  and  a 
wave  of  deficit  which  moves  to  the  left.  The  resulting 
height  field  has  a  height  gradient  in  geostrophic  balance 
with  the  velocity  in  the  region  of  the  jet  and  a  constant 
height  in  balance  with  the  zero  velocities  away  from  the 
jet.  Note  that,  unlike  initial  condition  1,  the 
balanced  state  requires  final  boundary  values  to  be 
different  from  their  initial  values.  Also  note  that  while 
the  Coriolis  effect  plays  little  role  in  the  evolution  of 
initial  condition  1,  it  is  essential  to  the  evolution  shown 
in  fig.  3.  If  f  =  0,  the  v  component  is  decoupled  from  the 
other  variables  in  equations  (1)  through  (3)  and  no 
adjustment  would  be  necessary. 

In  the  following  sections,  specific  lateral  boundary 
condition  formulations  will  be  described  and  the  results  of 
tests  using  these  conditions  in  the  one-dimensional  model 
will  be  presented.  All  boundary  conditions  were  tested  with 
both  initial  conditions.  In  the  figures  which  allow 
comparison  of  the  boundary  conditions,  the  height  field 
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ladjeled  "LD"  is  the  large  domain's  center  26  grid  points  and 
may  be  thought  of  as  the  "exact"  solution.  Each  boundary 
condition  formulation  will  be  denoted  with  a  three  or  four 
letter  descriptor  (such  as  FXD  for  fixed  conditions)  as  is 
used  in  the  figures. 


4.  Results  Using  Fixed  or  Extrapolated  Boundary  Conditions 


a)  Boundary  Condition  Formulations 

1)  Fixed  boundary  condition  (FXD).  This  boundary 
condition  consists  simply  of  specifying  boundary  values 
initially  and  holding  these  values  constant  throughout 
the  integration.  It  may  be  written 

where  <p  is  any  prognostic  variatble  i  rj  ,  u,  or  v  in 
the  present  model).  This  approach  eliminates  the 
potential  problem  of  geostrophically  adjusted  boundary 
errors,  provided  the  initial  boundary  values  are 
correctly  chosen.  This  boundary  condition  is  not 
mathematically  well-posed,  however,  in  that  it  results 
in  an  overspecification  of  the  problem.  It  also  is 
purely  reflective  of  waves  generated  in  the  domain's 
interior . 

2)  Combination  of  fixed  and  extrapolated  (FE2<).  In 
order  to  allow  the  boundary  to  be  more  "open",  while 


still  controlling  geostrophically  adjusted  errors,  some 
mesoscale  models  use  fixed  conditions  on  some  variaisles 
while  extrapolating  others  (Anthes  and  Warner,  1978; 
Nickerson,  1979).  The  typical  approach  is  to  fix  the 
pressure  and  thermodynamic  varicUbles  and  set  the 
horizontal  gradients  of  the  velocities  equal  to  zero 
(Anthes  and  Warner,  1978).  Here,  the  formulation 
follows  Nickerson  (1979)  with  the  height  specified  and 
the  velocities  extrapolated  only  when  the  flow  is 
outward.  That  is 


iiL  -  iv  .  « 

dx  "  dx  ^  for  u  outward  ,  (9) 

4^s4^'sO  foru  inward  . 

ot  ot 

3)  Zero-gradient  condition  (GRD) .  This  condition 
extrapolates  all  variables  on  the  boundary  by 
specifying  the  horizontal  gradient  to  be  zero.  This 
may  be  writtten 

=  0  ,  (10) 

where  <^  =  “Tj,  u,  or  v.  This  condition  is  implemented 
by  setting  the  boundary  value  equal  to  the  value  one 
grid  point  interior  to  the  boundary  at  each  timestep. 

4)  Sponge  condition  (SPG).  This  condition  follows  the 
formulation  of  Perkey  and  Kreitzberg  (1976).  The 


d</> 

dx 
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boundary  values  are  specified  and  do  not  change  in 
time.  However,  the  time  tendencies  at  the  gridpoints 
are  modified  such  that 

If  =  .  (11a) 

where 

0.0  for  I  =  boundary  grid  points 

0.4  for  I  =  boundary-1  grid  points 

W(I)  =  -0.7  for  I  =  boundary-2  grid  points  (11b) 

0.9  for  I  =  boundary-3  grid  points 

1.0  for  all  interior  grid  points 

This  type  of  boundary  condition  absorbs  waves  by 

reducing  their  phase  speed  to  zero  as  they  approach  the 

boundary  (Perkey  and  Kreitzberg,  1976).  As  the 

phasespeed  is  reduced,  the  wavelength  also  becomes 

smaller.  This  adds  energy  to  2  Ax  waves.  Therefore,  a 

filter  of  some  sort  is  required.  A  fourth-order 

diffusion  term  is  used  in  this  study  to  remove  this  2 Ax 

energy. 

b)  Comparison  of  the  Boundary  Condition  Tests 

Figures  4  and  5  show  the  height  field  in  comparative 
test  results  for  initial  condition  1  at  the  nondimensional 
times  t  =  36  and  t  =  72,  respectively.  The  height  field 
labeled  LD  is  the  center  26  grid  points  of  the  large  domain 
model  and  may  be  thought  of  as  the  solution  resulting  from 
truly  open  lateral  boundaries.  It  is  clear  that  all  of  the 


boundary  conditions  shown  here  have  serious  problems  with 
reflection  of  the  outward  propagating  waves. 

The  fixed  and  fixed/ extrapolated  (FXD  and  FEX) 
conditions  produced  nearly  identical  results.  This  suggests 
that  extrapolation  of  only  the  velocities  does  not  allow  the 
boundary  to  be  more  open  to  the  outgoing  wave. 

Extrapolation  of  all  the  variables  through  the  zero-gradient 
condition  (GRD)  is  not  adequate  either.  How<=ver, 
this  condition  appears  to  produce  errors  of  comparable 
magnitude  (though  of  different  phase).  The  sponge  condition 
while  the  best  of  this  set  of  conditions,  fails  to 
completely  absorb  the  wave.  This  may  be  a  result  of  the 
small  number  of  gridpoints  in  the  model  (which  is  typical  of 
the  number  of  grid  points  in  one  direction  in  a  3-D  model). 
It  may  also  be  a  result  of  the  relatively  rapid  phasespeed 
of  the  wave  in  this  model.  Tests  of  the  sponge  condition 
reported  by  Perkey  and  Kreitzberg  (1976)  which  showed  nearly 
completed  absorption  involved  internal  gravity  waves  which 
required  many  timesteps  to  traverse  the  width  of  the  sponge. 

A  similar  comparative  display  for  initial  condition  2 
at  times  t  =  36  and  t  =  72  is  shown  in  figures  6  and  7.  It 
is  impossible  for  any  of  the  boundary  conditions  which  hold 
the  boundary  value  fixed  (FXD,  FEX,  SPG)  to  achieve  the 
balanced  final  state.  This  leads  to  a  very  disturbed 
solution  as  the  initial  jet  in  the  v  component  continually 


forces  a  height  gradient  but  the  reflection  of  the  waves 
generated  interferes  with  the  height  field  adjustment.  The 
only  acceptable  balanced  final  state  for  these  conditions  is 
for  the  jet  in  v  to  be  diminished  to  zero.  Since  the 
geostrophlc  adjustment  on  this  scale  is  made  by  the  height 
field,  however,  v  is  only  diminished  by  the  small  viscosity 
in  the  model. 

The  zero-gradient  condition  (GRD)  which  is  not 
constrained  by  fixed  endpoint  values  appears  initially  to  be 
allowing  the  adjustment  to  proceed  with  only  small  errors. 

At  t  =  36,  the  height  gradient  in  the  jet  region  is 
apporximately  correct.  We  can  see  also  that  the  endpoints 
have  adjusted  in  the  right  sense,  although  they  are  further 
i  *om  the  equalibrium  height  than  they  should  be.  However, 
by  t  =  72,  it  is  clear  that  this  condition  is  not  handling 
the  adjustment  properly. 

Aqain,  the  sponge  condition  appears  to  yield  the  best 
results  despite  the  specified  values  on  the  endpoints.  The 
condition  tends  to  concentrate  the  imbalance  in  the  sponge 
region  near  the  boundaries  and  maintains  a  partial 
adjustment  in  the  interior.  While  not  ideal,  this  condition 
could  be  considered  adequate  in  many  situations,  especially 
if  the  initial  conditons  contained  only  small  imbalances. 


5.  Results  Using  Radiation  Type  Boundary  Conditions 


a)  Radiation  Condition  Formulations 

Radiation  type  lateral  boundary  conditions  are  based  on 
the  Sommerfeld  radiation  condition 

where  is  the  variable  of  interest  on  the  boundary  and  C 
is  the  phasespeed  of  the  outgoing  wave.  Variations  in  the 
radiation  condition  formulation  arise  from  the  various 
methods  that  can  be  used  to  specify  the  phasespeed,  C,  and 
to  finite  difference  the  equation.  For  the  shallow  water 
model,  all  waves  move  at  the  shallow  water  wavespeed,  so  C  = 
1.0  exactly.  In  the  general  case,  however,  the  wavespeed  is 
not  known  and  must  either  be  calculated  from  information 
available  in  the  interior,  or  specified  at  some  reasonable 
value.  Klemp  and  Lilly  (1978),  for  example,  specified  the 
phasespeed  as  the  expected  internal  gravity  wavespeed 
calculated  from  the  environmental  stability.  However,  Clark 
(1979)  showed  that  a  calculated  variable  phasespeed  gave 
less  reflection.  Most  studies  using  radiation  conditions 
have  used  a  variable  C.  We  will  consider  several  variable 
phasespeed  formulations. 

In  all  the  radiation  condition  formulations,  C  is 
required  to  satisfy  the  Courant-Fr iedrichs-Levy  iCFL) 
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d<^ 

dt 


condition.  In  addition,  if  C  is  calculated  to  be  inward,  it 
must  be  set  equal  to  zero  to  prevent  spurious  generation  of 
waves.  Thus,  C  must  satisfy 


0  for  C  inward 

C  for  0  <  C  <  Ax/At 

Ax/  At  for  C  >  Ax/At 


(13) 


Note  that  for  C  =  0,  (12)  reduces  to  d4>/dt  =  0  for  that 
timestep. 

1)  Orlanski-type  centered  on  previous  timestep  (ORP). 
This  is  the  original  formulation  by  Orlanski  (1976). 

It  calculates  C  by  inverting  (12)  and  evaluating 
d^/dt  and  d^/dx  at  the  previous  timestep  and  one  grid 
point  interior  to  the  boundary.  Then 


(14) 


where 


dt  •  2Ar[‘^b-i” • 


.  J,.,  ’  b-2  J  • 

In  the  above,  is  the  ith  grid  point  interior  to 

the  boundary  and  the  superscript  denotes  the  timestep. 
The  term  in  parenthesis  in  (16)  is  used  to  prevent 
tiraesplitting  on  the  boundary  point.  The  new  boundary 
value  is  then  given  by 
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^T»|  (i-CAt/Ax)  .T-I  .  2CA‘/Ax  .T 
■  D  ‘Pb  5 


{ 17) 


where  D  =  1  +  C  At /Ax. 

2)  Orlanski-type  centered  on  current  tiinestep  (ORC). 

As  suggested  by  Raymond  and  Kuo  (1984),  when  an 
explicit  time  scheme  is  used  in  a  model,  the 
calculation  of  C  can  be  centered  on  the  current 
timestep.  Then  the  calculation  of  C  is  identical  to 
(14)  through  (16)  with  r  replaced  by  t  +  1.  The 
boundary  value  calculation  is  still  identical  to  (17). 

3)  Upstream  technique  (UPS).  Bannon  (1979)  and  Miller 
and  Thorpe  (1981)  have  both  proposed  using  upstream 
differencing  to  evaluate  the  radiation  equation  (12). 
Inverting  the  equation  to  obtain  C  then  gives 


as  before,  but  now 


( 19) 


(20) 


The  boundary  point  is  then  given  by 


4)  Radiation  conditions  with  the  exact  phasespeed  (ORE 
and  UPE) .  For  the  shallow  water  model,  the  exact 
phasespeed  is  known  (C  =  1.0)  so  this  value  may  be  used 
in  the  radiation  condition  to  isolate  errors  resulting 
from  the  finite  difference  representation  alone.  When 
the  exact  phasespeed  is  used,  ORP  and  ORC  become 
identical  and  this  form  is  denoted  ORE.  The  upstream 
form,  UPS,  is  denoted  UPE  when  the  exact  phasespeed  is 
used.  This  approach  is  useful  in  the  present  study, 
but  is  not  applicable  in  a  more  general  model  which 
admits  waves  of  different  phasespeeds. 

b.  Comparison  of  the  Radiation  Condition  Tests 

Figures  8  and  9  show  a  comparison  of  the  radiation 
condition  formulations  for  initial  condition  1  at 
nondimens ional  times  t  =  36  and  t  =  72.  Again,  the  LD 
solution  may  be  thought  of  as  giving  an  exact  boundary 
value.  We  find  here  that  all  of  the  radiation  condition 
formulations  yield  superior  results  to  the  boundary 
conditions  used  to  produce  figures  4  and  5.  Differences 
between  the  individual  radiation  conditions  are  most 
apparent  at  t  =  72  (fig.  9). 

Clearly,  the  conditions  using  the  exact  phasespeed  are 
the  best  (ORE  and  UPE).  However,  even  these  result  in  a 
slight  reflection  of  wave  energy.  This  reflection  is  due  to 
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the  finite  differencing  of  (12)  and  shows  that  even  perfect 
calculation  of  C  in  a  more  general  model  will  result  in  a 
small  error  at  the  boundary.  This  result  is  also  noted  in 
the  formal  analysis  carried  out  by  Miller  and  Thorpe  (19B1). 

An  important  and  very  noticaQjle  point  in  the  three 
varicLble  phasespeed  formulation  results  (ORP,  ORC,  and  UPS) 
is  that  the  height  field  has  dropped  below  the  equilibrium 
level  of  17  =  1.0  at  t  =  72.  This  indicates  that  mass  has 
been  lost  through  the  open  boundaries  as  waves  exited.  This 
is  a  potentially  serious  problem  which  could  adversely 
affect  the  performance  of  a  more  complicated  atmospheric 
model.  This  is  clearly  related  to  the  variable  phasespeed 
since  the  exact  phasespeed  results  do  not  exhibit  a 
noticable  change  in  mass. 

The  flux  of  mass  through  the  boundaries  appears  to  be 
related  to  two  factors.  First,  the  phasespeeds  used  in  the 
radiation  condition  are  calculated  seperately  for  each 
variable.  Despite  the  strength  and  coherence  of  the  waves 
generated  by  initial  condition  1,  the  calculated  phasespeeds 
bear  little  resemblance  to  the  theoretical  ideal.  Fig.  10 
shows  the  phasespeeds  calculated  by  ORP  for  each  of  the 
three  variables  as  a  function  of  time.  As  can  be  seen,  they 
oscillate  wildy  between  zero  and  CFL  limit.  Further,  there 
is  little  relation  between  the  phasespeeds  calculated  for 
the  different  variables.  This  can  lead  to  boundary 


tendencies  which  are  not  consistent  with  continuity  of  mass 
and  a  resulting  mass  flux  through  the  boundary.  Plots  of 
the  phasespeeds  versus  time  for  ORC  and  UPS  show  similar 
characteristics. 

The  second  factor  which  appears  to  be  important  in  the 
flux  of  mass  through  the  boundary  is  the  nature  of  the 
staggered  grid.  The  boundary  gridpoints  on  which  the 
velocities  are  modified  by  the  boundary  condition  are  one 
half  grid  spacing  outward  from  the  boundary  points  for 
height.  Thus,  even  if  all  three  variaibles  use  the  same 
value  of  C,  this  value  is  being  used  at  different  spatial 
points  on  the  wave.  Thus,  unless  C  is  calculated  perfectly 
at  every  tiraestep,  different  parts  of  the  wave  will  be 
advected  at  different  speeds  through  the  boundary,  resulting 
in  a  net  mass  flux.  Results  obtained  using  averages  of  C, 
which  demonstrate  this  conclusion,  will  be  shown  in  the  next 
section. 

Figures  11  and  12  show  the  performance  of  the  radiation 
conditions  on  initial  condition  2  at  t  =  36  and  t  =  72. 

Since  the  boundary  values  are  able  to  evolve  naturally  in 
time,  all  of  the  radiation  condition  formulations  allow  a 
natural  adjustment  to  take  place.  Again,  ORE  and  UPE  give 
the  best  results,  but  all  of  the  formulations  give 
satisfactory  results.  There  was  little  evidence  of  a  change 
in  mass  in  the  experiments  using  initial  condition  2. 


c.  Results  Using  Averages  of  the  Phasespeeds 

As  can  be  seen  in  fig.  10,  while  the  individual 
phasespeeds  calculated  for  each  variable  oscillate  wildy,  an 
average  of  them  would  tend  to  be  closer  to  C  =  1.0.  Each  of 
the  radiation  condition  formulations  was  tested  with  the 
phasespeed  being  the  average  of  the  phasespeeds  calculated 
using  17,  u,  and  v  alone.  These  versions  of  ORP,  ORC,  and 
UPS  are  denoted  ORPA,  ORCA,  and  UPSA,  respectively.  The 
phasespeed  as  a  function  of  time  for  ORPA  using  initial 
condition  1  is  shown  in  fig.  10  as  the  heavy  solid  line. 

Note  that  this  is  not  exactly  the  same  as  the  average  of  the 
individual  phasespeeds  from  ORP  because  the  evolution  of  the 
flow  in  ORPA  is  modified.  As  anticipated,  the  phasespeed  in 
ORPA  (as  well  as  in  ORCA  and  UPSA  which  are  not  shown  but 
similar)  remained  closer  to  the  analytic  phasespeed,  C  = 

1.0. 

Figures  13  and  14  show  the  height  field  for  initial 
conditions  1  and  2  at  t  =  72  for  the  radiation  conditions 
using  the  average  phasespeed.  These  results  are  somewhat 
closer  to  the  results  obtained  with  the  analytic  phasespeed. 
It  is  evident  in  fig.  13  that  there  is  still  a  mass  flux 
through  the  boundary.  Now,  however,  the  mass  of  the  domain 
is  increasing  rather  than  decreasing. 

d.  The  Higher  Order  Upstream  Scheme  ( UPHO ) 

Miller  ar.i  Thorpe  (1981)  presented  a  radiation 


condition  formulation  using  upstream  differencing  which  has 
a  formal  accuracy  of  higher  order  than  ORP,  ORC,  or  UPS. 

The  model  results  using  this  condition  (denoted  UPHO)  with 
initial  condtions  1  and  2  are  shown  in  figures  13  and  14, 
respectively.  As  can  be  seen,  this  condition  does  not  allow 
the  waves  to  exit  the  domain  without  reflection.  Since  this 
condition  is  formally  more  accurate  than  the  previously 
discussed  radiation  conditions,  its  relatively  poor 
performance  is  somewhat  surprising.  Raymond  and  Kuo  (1984) 
tested  this  higher  order  scheme  in  a  two-dimensional  model 
and  also  found  it  to  give  unsatisfactory  results.  This 
section  will  outline  the  formulation  of  this  condition  and 
seek  to  explain  why  its  practical  accuracy  is  far  less  than 
its  formal  accuracy. 

Following  the  notation  of  Miller  and  Thorpe  (1981)  we 
can  let  r  =  C  At /Ax  be  a  nondimansional  phasespeed.  The 
CFL  and  ouflow  conditions  reduce  to  0  <  r  <  1  (for  r 
positive  on  outflow) .  Then  the  normal  upstream  calculation 
of  the  phasespeed  given  by  equaitons  (18)  through  (20)  may 
be  written 


L^b-l  ^b-l  JL^b- 
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(22) 


which  can  be  shown  to  be  of  second  order  accuracy  (Miller 
cmd  Thorpe,  1981)  Two  other  possible  calculations  of  r  are 


given  by 


and 


(23) 


(24) 


Miller  and  Thorpe  (1981)  then  show  that  the  proper 
combination  of  these  three,  given  by 

r  =  r|  r^  -  r^  (  25 ) 

achieves  a  formal  truncation  error  which  is  third  order. 
Thus,  this  combination  should  give  higher  accuracy.  To 
implement  this  condition,  we  find  r  using  equations  (22) 
through  (25),  then  let  C  =  r  Ax/At  and  use  equaiton  (21)  to 
update  the  boudary  value. 

Extensive  testing  with  the  one-dimensional  model  has 
confirmed  a  rather  simple  explanation  for  the  poor 
performance  of  this  condition.  The  calculations  of  r^  and  rj 
both  rely  on  an  upstream  difference  at  a  point  one  grid 
spacing  interior  to  the  boundary,  and  differ  only  in  the 
tiraestep  on  which  the  calculation  is  based.  For  an  existing 
wave  which  requires  several  timesteps  to  pass  through  the 
boundary,  both  r^  and  rj  would  be  expected  to  give  similar 
results.  The  calculation  of  r^,  however  depends  on  a  time 
difference  on  the  boundary  point  itself.  Therefore,  r^  is 
coupled  to  the  boundary  condition  itself  and  will  yield  a 
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nonzero  result  only  when  the  boundary  condition  allows  a 
change  in  the  boundary  value.  If  the  boundary  value  has  not 
changed  in  two  timesteps,  r^  =  0,  and  r  =  r^ -  r^ .  But 
r^  »  r^^ ,  so  their  difference  will  be  approximately  zero  and 
may  be  negative.  A  negative  r  indicates  inflow  (assuming  r 
positive  on  outflow) ,  so  r  is  set  to  zero  and  the  boundary 
value  is  not  modified.  This  feedback  maintains  rg  =  0  and 
results  in  essentially  fixed  conditions  with  the  resulting 
reflection.  In  the  simulations  made  using  UPHO,  brief 
periods  on  nonzero  r  interupted  longer  periods  when  the 
condition  set  r  to  zero.  Thus,  for  some  brief  periods  of 
the  integration,  waves  were  able  to  exit  the  domain. 

6.  Adjustment  to  Large  Scale  Imbalances 

The  results  for  initial  conditions  1  and  2  described 
above  involve  geostrophic  adjustment  of  imbalances  on  scales 
smaller  than  the  Rossby  radius  of  deformation.  In  these 
cases,  the  height  field  adjusts  to  the  winds.  The  model's 
response  to  imbalances  in  the  large  scale  ( "wavenumber 
zero")  fields  yields  results  of  importance  to  both  model 
initialization  and  lateral  boundary  condition  formulation. 
Even  though  the  model  domain  is  smaller  than  the  radius  of 
deformation,  a  constant  gradient  of  height  which  is  out  of 
balance  with  a  constant  wind  will  be  treated  by  the  model  as 
an  infinitely  large  disturbance.  In  this  case,  theory 


predicts  that  the  winds  should  adjust  to  the  pressure. 

One  experiment  involved  initializing  the  model  with  a 
constant  height  gradient  which  would  be  in  geostrophic 
balance  with  a  v  velocity  of  2  m/s.  The  initial  v  component 
was  set  to  zero,  however,  so  the  fields  were  out  of  balance. 
The  adjustment  process  in  this  case  is  fairly  clear  from  a 
consideration  of  equations  (1)  through  (3).  The  height 
gradient  accelerates  a  u  velocity  (from  (1))  which  in  turn 
accelerates  a  v  velocity  (from  (2)).  As  v  becomes  larger 
and  approaches  geostrophic  balance  with  the  height  gradient, 
the  acceleration  of  u  is  diminished  and  eventually  u  is 
decelerated  to  zero  leaving  only  the  geostrophically 
balanced  v  component  and  height  gradient.  Simulations  with 
the  one-dimensional  model  show  this  adjustment  process  to  be 
an  order  of  magnitude  slower  than  the  adjustment  of  a 
localized  disturbance.  By  the  nondimensional  time  t  =  240, 
the  V  component  had  been  accelerated  to  a  dimensional  value 
of  only  0.44  m/s.  Other  simulatior. s  using  a  variety  of 
large  scale  imbalanced  initial  conditions  have  confirmed 
this  slow  adjustment  process. 

This  result  indicates  that  it  is  very  important  for  the 
large  scale  gradients  to  be  very  close  to  balance  in  an 
atmospheric  model's  initialization.  The  local  imbalances 
due  to  topography  or  higher  resolution  data  present  in  the 
initialization  will  balance  quickly  provided  an  appropriate 
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lateral  boundary  condition  is  used.  Imbalances  in  the  large 
scale  gradients,  however,  will  adjust  in  a  time  scale 
comparable  to  the  inertial  period  ('-I?  hours). 

This  result  also  indicates  that  the  concern  of  Anthes 
and  Warner  (1978)  about  geostrophlcally  adjusted  boundary 
value  errors  is  not  warranted  in  meso-beta  models,  which 
typically  make  forecasts  for  periods  much  shorter  than  one 
day.  However,  it  may  still  be  important  for  larger  scale 
models  which  forecast  for  a  day  or  more.  The  adjustment 
process  is  sufficiently  slow  that  errors  produced  by  the 
boundary  conditions  will  require  several  hours  to 
signif icaintly  affect  the  mean  geostrophic  flow  in  the 
interior.  In  chosing  lateral  boundary  conditions  for 
meso-beta  models,  then,  emphasis  should  be  placed  on  the 
ability  of  the  condition  to  deal  with  waves  generated  in  the 
interior  of  the  domain  by  natural  processes  or  the 
geostrophic  adjustment  of  initial  local  imbalances. 

7.  Conclusions  and  Future  Activities 

This  report  has  discussed  the  test  results  lateral 
boundary  conditions  in  a  one-dimensional  shallow  water 
model.  The  results  show  conclusively  that  radiation 
conditions  allow  waves  to  exit  the  domain  without  serious 
reflection  and  allow  the  boundary  points  to  evolve  with  the 
geostrophic  adjustment  of  the  interior  so  that  a  balanced 


state  may  be  reached.  Of  the  non-radiation  type  conditions 
tested,  only  the  sponge  condition  of  Perkey  and  Kreitzberg 
(1976)  showed  signs  of  being  adequate,  though  still 
inferior.  One  of  the  primary  reasons  for  chosing  the  sponge 
condition  or  some  other  fixed  boundary  value  condition,  has 
been  the  concern  of  geostrophically  adjusted  large  scale 
gradient  errors  due  to  the  boundary  condition  (Anthes  and 
Warner,  1978).  This  has  been  shown  to  be  of  secondary 
importance  on  the  meso-beta  scale.  A  potential  problem  with 
radiation  conditions,  however,  is  the  flux  of  mass  through 
the  lateral  boundary  which  may  result  in  a  total  chamge  in 
mass  for  the  domain.  This  problem  is  reduced  if  the 
phasespeeds  are  calculated  accurately. 

The  analysis  of  lateral  boundary  conditions  using  a 
one-dimensional  model  does  not  replace  boundary  condition 
testing  in  a  full  three-dimensional  model.  It  does, 
however,  compliment  it.  The  one-dimensional  modeling  allows 
a  large  number  of  variations  of  boundary  conditions  to  be 
tested  economically.  Further,  since  the  one-dimensional 
model  equations  are  much  simpler  than  the  three-dimensional 
set,  physical  interpretation  of  the  results  is  much  simpler. 
The  price  of  this  simplification  is  that  some  types  of 
motion  (such  as  internal  gravity  waves)  may  be  filtered  out, 
so  the  easy  physical  interpretation  in  the  one-dimensional 
case  may  not  always  extend  to  the  three-dimensional  model 


results.  Thus,  It  Is  useful  to  use  both  approaches.  The 
one-dlmenslonal  model  can  eliminate  obviously  Inadequate 
boundary  conditions  and  test  minor  variations  in  the 
reasonable  ones,  and  the  three-dimensional  model  can  test 
the  most  promising  candidates  revealed  from  the 
one-dimensional  tests. 

Some  experiments  using  fixed/ extrapolated,  sponge,  and 
radiation  conditions  have  been  carried  out  with  the  full 
three-dimensional  model.  These  tests  have  basically 
supported  the  results  presented  in  this  report.  Future 
work  will  involve  further  testing  in  the  three-dimensional 
model.  Work  will  also  continue  with  the  one-dimensional 
model.  This  will  involve  the  inclusion  of  time  dependent 
data  on  the  lateral  boundaries  which  is  important  for  meso- 
beta  models  nested  in  larger  scale  models. 
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of  the  shallow  water  model  domain  showing' 


Figure  2.  Evolution  of  initial  condition  1  in  the  large  domain 
model.  Tick  marks  along  the  bottom  show  the  location  of 
grid  points,  and  the  dashed  lines  indicate  the  location  of 
endpoints  in  the  26  gridpoint  model. 


e  4.  Comparison  of  the  height  field  at  t  =  36  using  the 
indicated  boundary  conditions  in  the  26  gridpoint  model 
with  initial  condition  1. 


